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ABSTRACT 

Radio interferometry probes astrophysical signals through incomplete and noisy Fourier mea- 
surements. The theory of compressed sensing demonstrates that such measurements may ac- 
tually suffice for accurate reconstruction of sparse or compressible signals. We propose new 
generic imaging techniques based on convex optimization for global minimization problems 
defined in this context. The versatility of the framework notably allows introduction of spe- 
cific prior information on the signals, which offers the possibility of significant improvements 
of reconstruction relative to the standard local matching pursuit algorithm CLEAN used in 
radio astronomy. We illustrate the potential of the approach by studying reconstruction per- 
formances on simulations of two different kinds of signals observed with very generic inter- 
ferometric configurations. The first kind is an intensity field of compact astrophysical objects. 
The second kind is the imprint of cosmic strings in the temperature field of the cosmic mi- 
crowave background radiation, of particular interest for cosmology. 

Key words: techniques: interferometric, techniques: image processing, cosmology: cosmic 
microwave background 



1 INTRODUCTION 

Radio interferometry is a powerful technique for aperture syn- 
thesis in_astronomy, dating back to more than sixty years 
ago jRvle & Vonberd ImSt: iBlvt he 1957; iRvle et alj 1 19591 : 



. iRvle & HewishI Il960l : [Thompson et ah ,2004) . In a few words, 
' thanks to interferometric techniques, radio telescope arrays syn- 
thesize the aperture of a unique telescope of the same size as the 
maximum projected distance between two telescopes on the plane 
perpendicular to the pointing direction of the instrument. This al- 
. lows observations with otherwise inaccessible angular resolutions 
' and sensitivities in radio astronomy. The small portion of the celes- 
tial sphere accessible to the instrument around the pointing direc- 
tion tracked during observation defines the original real planar sig- 
nal or image I to be recovered. The fundamental Nyquist-Shannon 
theorem requires a signal to be sampled at a frequency of twice 
its bandwidth to be exactly known. The signal / may therefore be 
expressed as a vector x G containing the required number TV 
of sampled values. Radio-interferometric data are acquired in the 
Fourier plane. The number m of spatial frequencies probed may 
be much smaller than the number A'^ of discrete frequencies of the 
original band-limited signal, so that the Fourier coverage is incom- 
plete. Moreover the spatial frequencies probed are not uniformly 
sampled. The measurements are also obviously affected by noise. 
An ill-posed inverse problem is thus defined for reconstruction of 
the original image. 

Beyond the Nyquist-Shannon theorem, the emerging theory of 



com pressed sensing aims at merging data a cquisition and compres- 
sion dCandes et al.ll2006j|bl: lcmdes 2006; Donohdl2006l : lBaraniukl 
|2007|) . It notably relies on the idea that a large variety of signals in 
Nature are sparse or compressible. By definition, a signal is sparse 
in some basis if its expansion contains only a small number of non- 
zero coefficients. More generally it is compressible if its expansion 
only contains a small number of significant coefficients, i.e. if a 
large number of its coefficients bear a negligible value. Compressed 
sensing theory demonstrates that a much smaller number of linear 
measurements is required for accurate knowledge of such signals 
than is required for Nyquist-Shannon sampling. The sensing matrix 
must simply satisfy a so-called restricted isometry property. In par- 
ticular, a small number of random measurements in a sensing basis 
incoherent with the sparsity or compressibility basis will ensure this 
property with overwhelming probability, e.g. random Fourier mea- 
surements of a signal sparse in real or wavelet space. Consequently, 
if compressed sensing had been developed before the advent of ra- 
dio interferometry, one could probably not have thought of a much 
better design of measurements for sparse and compressible signals 
in an imaging perspective. 

In this work we present results showing that the theory of com- 
pressed sensing offers powerful image reconstruction techniques 
for radio-interferometric data. These techniques are based on global 
minimization problems, which are solved by convex optimization 
algorithms. We also emphasize on the versatility of the scheme rel- 
ative to the inclusion of specific prior information on the signal in 
the minimization problems. This versatility allows the definition of 
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image reconstruction techniques which are significantly more pow- 
erful than standard deconvolution algorithm called CLEAN used in 
the context of radio astronomy. 

In Section |2l we pose the inverse problem for image recon- 
struction from radio-interferometric data and discuss the standard 
image reconstruction techniques used in radio astronomy. In Sec- 
tion [3] we concisely describe the central results of the theory of 
compressed sensing regarding the definition of a sensing basis and 
the accurate reconstruction of sparse or compressible signals. In 
Section |4l we firstly comment on the exact compliance of radio 
interferometric measurements with compressed sensing. We then 
study the reconstruction performances of various compressed sens- 
ing imaging techniques relative to CLEAN on simulations of two 
kinds of signals of interest for astrophysics and cosmology. We fi- 
nally conclude in Section[5] 

Noti ce that a first app lication of compressed sensing in as- 
tronomy teobin et al.ll2008l ) was very recently proposed for non- 
destructive data compression on board the future Herschel space 
observatorjQ. The versatility of the compressed sensing frame- 
work to account for specific prior information on signals was 
already pointed out in that context. Moreover, the generic po- 
tential of compressed sensing for interferometry was pointed in 
the signal processing community since the t ime when the the- 
ory emerged yPonoho 2006; Cand es et al.ll2006 b: Mary & Michel 
l2007l : iLevanda & Leshenj 1200^ It was also very recently ac- 
knowledged in radio astronomy 1comwelll2008h . The present work 
nonetheless represents the first application of compressed sens- 
ing for the definition of new imaging techniques in radio in- 
terferometry. A huge amount of work may be envisaged along 
these lines, notably for the transfer of the proposed techniques 
to optical and infrared interferometry. The extension of these 
techniques from the plane to the sphere will also be essential, 
notably with regard to forthcoming radio interferometer s with 
wide fields of view on t he celestial sphere (Comwell et ai] |2008l : 
iMcEwen & Scaifell2008h. such as the future Square Kilometer Ar- 
ray (SKAf I jCarilli & Rawlingsll2004h . 



2 RADIO INTERFEROMETRY 

In this section, we recall the van Cittert-Zernike theorem on the ba- 
sis of which we formulate the inverse problem posed for image re- 
construction from radio-interferometric data. We also describe and 
discuss the standard image reconstruction techniques used in ra- 
dio astronomy, namely a local matching pursuit algorithm called 
CLEAN and a global optimization algorithm called the maximum 
entropy method (MEM). 

2.1 van Cittert-Zernike theorem 

In a tracking configuration, all radio telescopes of an interferomet- 
ric array point in the same direction. The field of view observed on 
the celestial sphere S^ is limited by a so-called illumination func- 
tion A{u]), depending on the angular position a; G S^. The size of 
its angular support is essential ly inversely proportion al to the size 
of the dishes of the telescopes ( [Thompson et al.|[2004t) . At each in- 
stant of observation, each telescope pair identified by an index b 
measures a complex visibility i/;, £ C. This visibility is defined as 

^ http://herschel.esac.esa.int/ 
^ http://www.skatelescope.org/ 



the correlation between incoming electric fields E at the positions 
of the two telescopes in the three-dimensional space, 6i, &2 G R'': 

2/6= (&2,t) ^ . (1) 

In this relation, t denotes the time variable and the brackets (•)At 
denote an average over a time At long relative to the period of the 
radio wave detected. 

We consider a monochromatic signal with a wavelength of 
emission A, and made up of incoherent sources. We also consider a 
standard interferometer with an illumination function whose angu- 
lar support is small enough so that the field of view may be identi- 
fied to a planar patch of the celestial sphere: P C K^. The signal 
and the illumination function thus respectively appear as functions 
I(p) and A(p) of the angular variable seen as a two-dimensional 
vector p £ with an origin at the pointing direction of the array. 
The vector Bb = &2 — bi G R^ defining the relative position be- 
tween the two telescopes is called the baseline, and its projection on 
the plane perpendicular to the pointing direction of the instrument 
may be denoted as £ R^. One also makes the additional as- 
sumption that the maximum projection of the base lines in the point- 
ing direction itself is small jCornwell et al.ll20o3) . In this context, 
the so-called van Cittert-Zernike theorem states that the visibility 
measured identifies with the two-dimensional Fourier transform of 
the image multiplied by the illumination function AI at the single 
spatial frequency 



i.e. 

yb = AI (ut ) , (3) 
with 

AI{u)= f A{p)I{p)e-^'^^^-^d'p, (4) 

for any two-dimensional vector u £ R^. Interferometric arrays thus 
probe signals at a resolution equivalent to that of a single tele- 
scope with a size R essentially equivalent to the maximum pro- 
jected baseline on the plane perpendicular to the pointing direction: 
R ~ maxj, B^. This expresses the essence of aperture synthesis 
dThompson et al.ll2004l) . 

2.2 Interferometric inverse problem 

In the course of an observation, the projected baselines on the plane 
perpendicular to the pointing direction change thanks to the Earth's 
rotation and run over an ellipse in the Fourier plane of the origi- 
nal image, whose parameters are linked to the parameters of ob- 
servation. The total number m/2 of spatial frequencies probed by 
all pairs of telescopes of the array during the observation provides 
some Fourier coverage characterizing the interferometer. Any in- 
terferometer is thus simply identified by a binary mask in Fourier 
equal to 1 for each spatial frequency probed and otherwise. The 
visibilities measured may be denoted as a vector of m/2 complex 
Fourier coefficients y G _ _ AI{ub)}i^b!im/2, pos- 

sibly affected by complex noise values n G C™'^^ = {nt — 
"-('"6)}i^b^m/2 of astrophysical or instrumental origin. Consid- 
ering that the signal / and the illumination function A are real, 
a symmetry AI{—ut) = AI (ui,) also holds so that indepen- 
dent measurements may all be localized in one half of the Fourier 
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plane. The binary mask in Fourier identifying tlie interferome- 
ter is defined in this half of the plane and rendered symmetric 
around the origin so that it also corresponds to the Fourier trans- 
form of a real function. In this context, the measured visibilities 
may equivalently be denoted as a vector of m real Fourier coeffi- 
cients y £ R"* = {yrji^rsj m consisting of the real and imaginary 
parts of the complex measures, possibly affected by real noise val- 
ues n G = {jlrjl^rsSm. 

The original signal I(p) and the illumination function A(p) 
can be approximated by band-limited functions restricted to the fi- 
nite field of view precisely set by the illumination function: p £ P. 
In this context, we notice that they are identified by their Nyquist- 
Shannon sampling on a discrete uniform grid of A'^ = N^^^ x N^^'^ 
in real space with 1 ^ i ^ N . The sampled 



= {Xi = I{Pi)}l^i^N 

= {ai = 



points Pi £ 

signal may thus be denoted as a:: G 
while the illumination function is denoted as a G 
A{pi)}\^i^N , and the sampled product reads as x G — {xi = 
AI{pi)}i^i^M ■ Because of the assumed finite field of view, the 
functions may equivalently be described by their complex Fourier 
coefficients on a discrete uniform grid of N = N^^'^ x A''^''^ spa- 
tial frequencies Ui with 1 ^ i ^ A''. This grid is symmetric around 
the origin and limited at the maximum frequency defining the band 
limit. In particular for the Fourier coefficients of the product AI 
one has: i G = {Si = AI{ui)} i<i<N- The functions being 
real, one again has the symmetry AI{—Ui) — AI (ui) so that the 
signal is described by exactly N/2 complex Fourier coefficients in 
one half of the Fourier plane, or equivalently A*' real Fourier coef- 
ficients consisting of the real and imaginary parts of these complex 
coefficients. In the following we only use this decomposition with 
real coefficients in one half of the Fourier plane. 

However the frequencies Ub probed defined by (|2} for 1 ^ 
6 ^ m/2 are continuous and do not generally belong to the set of 
discrete frequencies Ui for 1 ^ i ^ A^. Reconstruction schemes in 
general perform a preliminary gridding operation on the visibilities 
yr with 1 ^ r ^ m so that the inverse problem may be refor- 
mulated in a pure discrete setting, i.e. between the discrete Fourier 
and real planes ( Thompson et al. 2004). The essential reason for the 
gridding resides in the subsequent use of the standard fast Fourier 
transform (FFTfl. For the sake of the considerations that follow we 
assume that the frequencies probed ut belong to the discrete grid 
of points Ui so that no artifact due to the gridding is introduced. In 
this discrete setting the Fourier coverage is unavoidably incomplete 
in the sense that the number of real constraints m is always smaller 
than the number of unknowns N: m < N. An ill-posed inverse 
problem is thus defined for the reconstruction of the signal x from 
the measured visibilities y as: 



y — ^^.x + n, 



(5) 



for a given noise n, and with a sensing matrix for radio inter- 
ferometry of the form 



MFD. 



(6) 



In this relation, the matrix D G R^'''^ = {Da, = 
<iiSii'}i!ii,i' s^N is the diagonal matrix implementing the illumina- 
tion function, and the matrix F G R'^^^ = {-Fii' }i^i,i':g]v im- 
plements the discrete Fourier transform providing the real Fourier 
coefficients in one half of the Fourier plane. The matrix M G 



^ Notice that fast algorithms have been developed to c ompute a Fourie r 
transform on non-equispaced spatial frequencies (NFFT) jPotts et alj2008l) . 
This could in principle allow one to avoid an explicit gridding operation. 



jjmxAf _ {j\,fri}is;,<m;isji<]v is the fectangular binary matrix 
implementing the mask characterizing the interferometer in one 
half of the Fourier plane. It contains only one non-zero value on 
each line, at the index of one of the two real Fourier coefficients 
corresponding to each of the spatial frequencies probed. 

We restrict our considerations to independent Gaussian noise 
with variance af. = a^{yr). From a statistical point of view, the 
likelihood £ associated with a candidate reconstruction x* of the 
signal x is defined as the probability of the data y given the model 
X*, or equivalently the probability of the noise residual n* — y — 
Under the Gaussian noise assumption it reads as 



C {y\x*) oc exp 



-X^ (x*;$;,t/) 



with the corresponding negative logarithm 



(a;*;i'„,y) = E 



«) 



(7) 



(8) 



following a chi-square distribution with m degrees of freedom. The 
defines a noise level estimator. The level of residual noise n* 
should be reduced by finding x* minimizing this x^, which corre- 
sponds to maximize the likelihood C. Typically, the measurement 
constraint on the reconstruction may be defined as a bound 



(9) 



with corresponding to some (lOOa)'^ percentile of the chi- 
square distribution, i.e. p{x^ ^ e'^) = a for some a 1. For 
a solution with a = there is a probability a that pure noise 
gives a residual smaller than or equal to the observed residual n* , 
and a probability 1 — a that noise gives a larger residual. Too small 
an a would thus induce possible noise over-fitting, i.e. inclusion of 
part of the noise in the reconstruction. These considerations might 
of course be generalized to other kinds of noise distributions. 

The inverse problem being ill-posed, many signals may for- 
mally satisfy measurement constraints such as l[9]l. In general, the 
problem may only find a unique solution a::*, as close as possible to 
the true signal x, through a regularization scheme which should 
encompass enough prior information on the original signal. All 
possible image reconstruction algorithms will essentially be dis- 
tinguished through the kind of regularization considered. 

2.3 Standard imaging techniques 

The general inverse problem ^ is to be considered if one wishes to 
undo the multiplication by the illumination function and to recover 
the original signal x on the given field of view. In practice, the re- 
construction is usually considered for the original image / already 
multiplied by the illumination function A, whose sampled values 
are x — Dx G R^ = {xi — aiXi}i^i^isr ■ In this setting the 
inverse problem reads as 



y = ^^^x + n, 



(10) 



with a sensing matrix strictly implementing a convolution: 

^^,=MF. (11) 

Firstly, the most standard and otherwise already very effec- 
tive image reconstruction algorithm from visibility measurements 
is called CLEAN. It approaches the image reconstruction in te rms 
of the corresponding deconvolution problem in real space (H ogbomI 
Il974l : ISchwarz|[T97a: iThompson et aLll2004l) . In standard vocabu- 
lary, the inverse transform of the Fourier measurements with all 
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non-observed visibilities set to zero is called the dirty image. Its 
sampled values a;''*' G = {x^^^}i^i^N are simply obtained 
by application of the adjoint sensing matrix to the observed visibil- 



on the signal (lAblesll974lGull & Daniellll999l : ICornwell & EvansI 



I1985| : IGu11 & Skillinelll999h . For positive signals, the relative en- 
tropy function between a sampled signal x £ = {xiji^i^jv 



ities: S*'*' = ^J.y. The inverse transform of the binary mask identi- and a model 2: G R^ = {zi}i^is;;v takes the simple form 



fying the interferometer is called the dirty beam. Its sampled values 
d e R'^ = {diji^i^jv follow from the application of the adjoint 
sensing matrix to a vector of unit values Im G R™: d = ^\ Im- 
The inverse transform of the noise n with all non-observed visibili- 
ties set to zero defines an alternative expression of the noise in real 
space. Again its sampled values rz''*' G R^ = {n'''^ ji^i^jv are 
simply obtained by application of the adjoint sensing matrix to the 
noise realization: n''*' = ^^.n. The inverse problem ilOi can thus 
be rephrased by expressing the dirty image as the convolution of 
the original image with the dirty beam, plus the noise: 

^d*x + n^'^K 



(12) 



CLEAN is a non-linear deconvolution method relying on this re- 
lation and working by local iterative beam removal. At each iter- 
ation, the point in real space is identified where a residual image, 
initialized to the dirty image, takes its maximum absolute value. 
The beam is removed at that point with the correct amplitude to 
produce the residual image for the next iteration. Simultaneously 
the maximum absolute value observed renormalized by the central 
value of the beam is added at the same point in the approximation 
image, initialized to a null image. This procedure assumes that the 
original signal is a sum of Dirac spikes. A sparsity or compressibil- 
ity prior on the original signal in real space is implicitly introduced 
so that its energy is concentrated at specific locations. On the con- 
trary, the Gaussian noise should be distributed everywhere on the 
image and should not significantly affect the selection of points in 
the iterations. This underlying sparsity hypothesis serves as a regu- 
larization of the inverse problem. 

A loop gain factor 7 is generally introduced in the procedure 
which defines the fraction of the beam considered at each itera- 
tion. Values 7 around a few tenths are usually used which allow for 
a more cautious consideration of the sidelobes of the dirty beam. 
The overall procedure is greatly enhanced by this simple improve- 
ment, albeit at high computational cost. In a statistical sense, the 
stopping criterion for the iteration procedure should be set in terms 
of relation l|9ll. However, the procedure is known to be slow and the 
algorithm is often stopped after an arbitrary number of iterations. 

Various weighting schemes can be applied to the binary mask 
in Fourier. Natural weighting simply corresponds to replace the unit 
values by the inverse variance of the noise affecting the correspond- 
ing visibility measurement. This corresponds to a standard matched 
filtering operation allowing the maximization of the signal-to-noise 
ratio of the dirty image before deconvolution. So-called uniform 
and robust weightings can notably be used to correct for the non- 
uniformity of the Fourier coverage associated with the measured 
visibilities and to reduce the sidelobes of the dirty beam in real 
space. Multi-sca le versions of this method were also developed 
dComwelfcOOSh . 

CLEAN and multi-scale versions may actually be formu- 
lated in terms of the well-known matching pursuit (MP) procedure 
iMallat & Zhang 1993; Mallat 1998). The corresponding MP algo- 
rithm simply uses a circulant dictionary for which the projection 
on atoms corresponds to the convolution with the dirty beam. The 
loop gain factor may also be trivially introduced in this context. 

Secondly, another approach to the reconstruction of images 
from visibility measurements is MEM. In contrast to CLEAN, 
MEM solves a global optimization problem in which the inverse 
problem ilOl is regularized by the introduction of an entropic prior 



5* (x, 2) = — Xi In — . 



(13) 



This function is always negative and takes its maximum null value 
when X — z.In the absence of a precise knowledge of the signal x, 
z is set to a vector of constant values. In such a case, maximizing 
the entropy prior promotes smoothness of the reconstructed image. 

The MEM problem is the unconstrained optimization problem 
defined as the minimization of a functional corresponding to the 
sum of the relative entropy 5* and the '■ 



mm 



(14) 



for some suitably chosen regularization parameter r > 0. In gen- 
eral, the minimization thus requires a trade-off between mini- 
mization, and relative entropy maximization. 

Notice that the definition il3\ may easily be generalized for 
non-positive signals. A multi-scale version of MEM was also de- 
fined. It considers that the original image may have an efficient 
representation in terms of its decomposition in a wavelet basis. The 
entrop y is then defined direc tly on the wavelet coefficients of the 
signal ( Maisinge r et al 1 1999 *). 

For completeness we finally quote the WIPE reconstruction 
procedure which also solves a global minimization problem, but in 
which the inverse problem JlOb is regularized by the introduction 
of a smoothness prior on the part of the signal whose Fourier sup- 
port corresponds to the non-probed spatial frequencies. This corre- 
sponds to minimize the y* ^ after assigning a null valu e to all initially 
non-observed visibilities jLannes et al ]|l994ll996l) . 

In conclusion, CLEAN is a local iterative deconvolution tech- 
nique, while MEM and WIPE are reconstruction techniques based 
on global minimization problems. All three approaches are flexi- 
ble enough to consider various bases (Dirac, wavelet, etc.) where a 
majority of natural signals can have a sparse or compressible repre- 
sentation. CLEAN also implicitly assumes the sparsity of the signal 
in the reconstruction procedure. But none of these methods explic- 
itly imposes the sparsity or compressibility prior on the reconstruc- 
tion. This precise gap is notably bridged by the imaging techniques 
defined in the framework of the compressed sensing theory. 



3 COMPRESSED SENSING 

In this section we define the general framework of the theory 
of compressed sensing and quote its essential impact beyond the 
Nyquist-Shannon sampling theorem. We then describe the re- 
stricted isometry property that the sensing basis needs to satisfy 
so that sparse and compressible signals may be accurately recov- 
ered through a global optimization problem. We finally discuss the 
idea that incoherence of the sensing and sparsity or compressibil- 
ity bases as well as randomness of the measurements are the key 
properties to ensure this restricted isometry. 



3.1 Beyond Nyquist-Shannon 

In the framework of compressed sensing the signals probed are 
firstly assumed to be sparse or compressible in some basis. Tech- 
nically, we consider a real signal identified by its Nyquist-Shannon 
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sampling as s G = {a:^j}i^i^jv. A real basis * G R^^^ = 
{"J/iiuji^isj iva^jii^T is defined, which may be either orthog onal, 
with r = iV, or redundant, with T > N jRauhut et al.ll2008l) . The 
= {a™}i^iD^T of the signal defined by 



decomposition a G 



(15) 



is sparse or compressible in the sense that it only contains a small 
number K <^ N of non-zero or significant coefficients respec- 
tively. The signal is then assumed to be probed by m real linear 
measurements y G R™ = {i/r}i^r^m in some real sensing basis 
ir^r^m;i^i^N and possibly affected by inde- 
pendent and identically distributed noise n G R™ = {nrjisjrsSm: 

y = ea + n withe = G R™""^. (16) 



This number m of constraints is typically assumed to be smaller 
than the dimension A'^ of the vector defining the signal, so that the 
inverse problem ( I16t is ill-posed. 

In this context, the theory of compressed sensing defines the 
explicit restricted isometry property (RIP) that the matrix Q should 
satisfy in ord er to allow an accurat e recovery of spa rse or compress- 
ible signals dCandes et al.] 12006 j|bl ; ICandesl |2006l). In that regard, 
the theory offers multiple ways to design suitable sensing matri- 
ces $ from properties of incoherence with and randomness of 
the measurements. It shows in particular that a small number of 
measurements is required relative to a naive Nyquist-Shannon sam- 
pling: m <^ N. The framework also defines a global minimiza- 
tion problem for the signal recovery called Basis Pursuit (BP). This 
problem regularizes the originally ill-posed inverse problem by an 
explicit sparsity or compressibility prior on the signal. The corre- 
sponding solution may be obtained through convex optimization. 
Alternative global minimization problems may also be designed. 



3.2 Restricted isometry and Basis Pursuit 

Let us primarily recall that the £p norm of a real vector u G 

C*^ — {ui}i^i^Q is defined for any p G R+ as ||it||p = 



|p\i/j 



, where stands for the absolute value of the 



component ui . The well-known £2 norm is to the square-root of the 
sum of the absolute values squared of the vector components. 

By definition the matrix O satisfies a RIP of order K if there 
exists a constant Sk < 1 such that 



(1 - Sk) ||qk||2 < ||0aif II2 ^ (1- 5b 



\aK\\l, 



(17) 



for all vectors qa' containing at maximum K non-zero coefficients. 

The £1 norm of the vector a G R'^ — {amji^m^T is simply 
defined as the sum of the absolute values of the vector components: 



Q 1 



(18) 



From a Bayesian point of view, this £1 norm may be seen as the 
negative logarithm of a Laplacian prior distribution on each inde- 
pendent component of a. For comparison the square of the £2 norm 
may be seen as the negative logarithm of a Gaussian prior distribu- 
tion. It is well-known that a Laplacian distribution is highly peaked 
and bears heavy tails, relative to a Gaussian distribution. This cor- 
responds to say that the signal is defined by only a small number of 
significant coefficients, much smaller than a Gaussian signal would 
be. In other words the representation a of the signal x in the spar- 
sity or compressibility basis ^ is indeed sparse or compressible if 
it follows such a prior. Finding the a' that best corresponds to this 
prior requires to maximize its Laplacian probability distribution, or 



equivalently to minimize the £1 norm. Notice that th is conclusion 
also follows from a pur e geometrical argument in R^ jCandes et al.l 
l2006bl : lBaraniugi2007l) . 

A constrained optimization problem explicitly regularized by 
a £1 sparsity prior can be defined. This so-called Basis Pursuit de- 
noise (BPe) problem is the minimization of the £1 norm of a' under 
a constraint on the £2 norm of the residual noise: 



min a 1 subject to 



\y — Qa'\\2 ^ 



(19) 



Let us recall that the noise was assumed to be identically dis- 
tributed. Consequently, considering Gaussian noise, the £2 norm 
term in the BPe problem is identical to the condition {9}, for 
corresponding to some suitable percentile of the distribution 
with m degrees of freedom governing the noise level estimator. 

This BPe problem is solved by application of non-linear and iter- 

" 71 

ative convex optimization algorithms (Combettes & Pesquet 200g; 

van den Berg & Friedlander 2008). In the absence of noise, the BPe 
problem is simply called Basis Pursuit (BP). If the solution of the 
BPj problem is denoted a* then the corresponding synthesis-based 
signal reconstruction reads, from dlSt . as x* — 'ia* . 

Compressed sensing shows that if the matrix Q satisfies a RIP 
of ord er 2K with some suitable constant S2K < \/2 — 1 cCandfesI 
l2008l) . then the solution 2:* of the BPj problem provides an accurate 
reconstruction of a signal x that is sparse or compressible with K 
significant coefficients. The reconstruction may be said to be opti- 
mal in that exactly sparse signals are recovered exactly through BP 
in the absence of noise: x* = x. Moreover strong stability results 
exist for compressible signals in the presence of noise. In that case, 
the £2 norm of the difference between the representation a of the 
signal in the sparsity or compressibility basis and its reconstruction 
a* is bounded by the sum of two terms. The first term is due to 
the noise and is proportional to e. The second term is due to the 
non-exact sparsity of a compressible signal and is proportional to 
the £1 norm of the difference between a and the approximation uk 
defined by retaining only its K largest components and sending all 
other values to zero. In this context, one has 



K 



(20) 



for two known constants Ci,k and C2,k depending on 52k- For 
inst ance, when 52k = 0-2, we have C i.k ~ 8.5 and C2.K ~ 
4.2 jCandes et alj|2006bl : ICandfesll2008l) . In an orthonormal basis 
this relation represents an explicit bound on the £2 norm of the 
difference between the signal x itself and its reconstruction x* as 
— a;*jl2 ~ j|a — Q*j|2. Moreover a;^ — ^'qa' then represents 
the best sparse approximation of x with K terms, in the sense that 
1 1 x — a; A- 1 1 2 is minimum. 

The constrained BPe problem may also be rephrased in terms 
of an unconstrained minimization problem for a functional defined 
as the sum of the £1 norm of a' and the £2 norm of the residual 
noise: 



mm 



y 



fc)a 1 12 + t| |a 



(21) 



for some suitably chosen regularization parameter r > 0. For 
each value of e, there exists a value r such that the solutions of 
the constrained and unconstrained £1 s parsity problems are identi- 
cal ( I van den Berg & FriedlandedboOSh . From a Bayesian point of 
view, this minimization is then equivalent to maximum a posteriori 
(MAP) estimation for a signal with Laplacian prior distribution in 
the sparsity or compressibility basis, in the presence of Gaussian 
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Finally, alternative minimization problems may be defined for 
the recovery. Firstly, a Ip norm with < p ^ 1 may for example 
be substituted for the li norm in the definition of the minimization 
problem. From a Bayesian point of view, the Ip norm to the power p 
may be seen as the negative logarithm of a prior distribution iden- 
tified as a generalized Gaussian distribution (GGD). Such distri- 
butions are even more highly peaked and bear heavier tails than a 
Laplacian distribution and thus promote stronger compressibility of 
the signals. Theoretical results hold for such Ip norm minimization 
problems when a RIP is satisfied ( F oucart & Lail2008 ). Such prob- 
lems are non-convex but can be solved iteratively by convex op- 
timi zation algorithms performing re-weighted ii norm minimiza- 
tion teandes etal."2008'; 'Pavi es & Gribonvalll2008l : iFoucart & Lail 
I2OO8 : Chartrand & Yin 2007). Secondly, a TV norm may also be 
substituted for the l\ norm in the definition of the minimization 
problem for signals with sparse or compressible gradients. The TV 
norm of a sign al is simply defined as the l\ norm of the magnitude 
of its gradient jRudin et al.lll992 ^. A theoretical result of exact re- 
construction holds for such TV norm minimization problems in the 
case of Fourier measurements of si gnals with exact ly sparse gradi- 
ents in the absence of noise ( iCandes et al. I l2006al) . But no proof 
of stability relative to noise and non-exact sparsity exists at the 
moment. Such minimization is also accessible through an itera tive 
schem e from convex optimization algorithms ( Candes & Romberg 
l2005h . 

This flexibility in the definition of the optimization problem is 
a first important manifestation of the versatility of the compressed 
sensing theory, and of the convex optimization scheme. It opens the 
door to the definition a whole variety of powerful image reconstruc- 
tion techniques that may take advantage of some available specific 
prior information on the signal under scrutiny beyond its generic 
sparsity or compressibility. 



3.3 Incoherence and randomness 

The issue of the design of the sensing matrix <I> ensuring the RIP 
for 6 = is of course fundamental. One can actually show that 
incoherence of $ with the sparsity or compressibility basis 4' and 
randomness of the measurements will ensure that the RIP is sat- 
isfied with overwhelming probability, provided that the number of 
measurements is large enough relati ve to the sparsity K considered 



fl = V N max I {fe\lpe' ) I 



The RIP is then satisfied if 



Cm 
^2 In* iV ' 



(23) 



(24) 



for some constant C' . As the incoherence is maximum between the 
Fourier and real spaces with /i = 1, the lowest number of measure- 
ments would be required for a signal that is sparse in real space. 
Notice that a factor In N instead of In* A'^ in condition ( 1241 ) was not 
proven but conjectured, suggesting that a lower number of mea- 
sure ments would still en sure the RIP. In that regard, empirical re- 
sults jLustig et suggest that ratios mjK between 3 and 5 
already ensure a reconstruction quality through BPe that is equiva- 
lent to the quality ensured by l |20t . 

Let us also emphasize that the TV norm minimization is often 
used from Fourier measurements of signals with sparse or com- 
pressible gradients. As already stated no stability result such as 
( 1201 ) was proven for the reconstruction provided by this minimiza- 
tion scheme. Empirical results suggest however that TV norm min- 
imization provides the same quality of recon struction as BPg for 
the same typical ratios m / K between 3 and 5 iCandes & RomberS 
l2005l : lLustigetalj2007l) . 



4 APPLICATIONS 

In this section, we firstly comment on the exact compliance of radio 
interferometric measurements with compressed sensing. We then 
consider simulations of two kinds of signals for reconstruction from 
visibility measurements: an intensity field of compact astrophysical 
objects and a signal induced by cosmic strings in the temperature 
field of the cosmic microwave background (CMB) radiation. Rely- 
ing on the versatility of the convex optimization scheme, enhanced 
minimization problems are defined in the compressed sensing per- 
spective through the introduction of specific prior information on 
the signals. The reconstruction performance is studied in compar- 
ison both with the standard BPe reconstructions in the absence of 
specific priors and with the CLEAN reconstruction. 



measurements is large enough relative to trie sparsity K considered a ^ ^ ^ c <-j j- 

1 II — - — H I — " II n ^ ^ 4.1 Interferometric measurements and compressed sensing 

jCandes et al.ll2006bl : ICandesll200a) . In this context, the variety of 



approaches to design suitable sensing matrices is a second form of 
the versatility of the compressed sensing framework. 

As a first example, the measurements may be drawn from a 
Gaussian matrix $ with purely random real entries, in which case 
the RIP is satisfied if 



Cm 



ln(iV/m) 



(22) 



for some constant C. The most recent result provides a value C ~ 
0.5, hence showing that the required redunda ncy of measurements 
mjK is very small jPonoho & TanneiJl2009l) . 

As a second example of interest for radio interferometry, 
the measurements may arise from a uniform random selection of 
Fourier frequencies. In this case, the precise condition for the RIP 
depends on the degree of incoherence between the Fourier basis 
and the sparsity or compressibility basis. If the unit-normed basis 
vectors corresponding to the lines of F and the columns of '3/ are 
denoted {/eji^e^jv and {i/'e' }i^e'^T> the mutual coherence /i of 
the bases may be defined as their maximum scalar product: 



In the context of compressed sensing, the sensing matrix needs to 
satisfy the RIP. If Fourier measurements are considered, this re- 
quirement may be reached through a uniform random selection of 
a low number of Fourier frequencies. In the context of radio in- 
terferometry, realistic visibility distributions are deterministic, i.e. 
non-random, superpositions of elliptical distributions in the Fourier 
plane of the image to reconstruct. However, the structure of the 
Fourier sampling is extremely dependent on the specific configu- 
ration of the radio telescope array under consideration. Visibilities 
from various interferometers may be combined, as well as visibil- 
ities from the same interferometer with different pointing direc- 
tions in the mosaicking technique (Thompson et al. 2004). From 
this point of view the realistic visibility distributions themselves 
are rather flexible. Moreover, the standard uniform weighting of the 
visibilities may be used to provide uniformity of the effective mea- 
surement density in the Fourier plane. Correctly studied realistic 
distributions might thus not be so far from complying exactly with 
the compressed sensing requirements. Finally, it was recently sug- 
gested that specific deterministic distributions of a low number of 
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linear measurements might in fact allow ac curate signal reconstru c- 
tion in the context of compressed sensing jMatei & Meve3l2008l) . 

Nonetheless, modifications of radio interferometric measure- 
ments might be conceived in order to comply exactly with stan- 
dard compressed sensing results. To this end, one might want to 
introduce randomness in the visibility distribution. Formally, ran- 
dom repositioning of the telescopes during observation or random 
integration times for the definition of individual visibilities could 
provide important advances in that direction. Also notice that com- 
pressed sensing does not require that measurements be identified to 
Fourier coefficients of the signal. The versatility of the framework 
relative to the design of suitable sensing matrices might actually 
be used to define generalized radio interferometric measurements, 
beyond standard visibilities, ensuring that the RIP is explicitly sat- 
isfied. In this perspective, direct modifications of the acquisition 
process through a scheme similar to spread spe ctrum techniques 
jNaini et al.ll2009l) or coded aperture techniques jMarcia & Willetll 
l2008l) could also provide important advances. 

In the following applications we simply consider standard vis- 
ibility measurements. We assume generic interferometric configu- 
rations characterized by uniform random selections of visibilities. 

4.2 Experimental set up 

We consider two kinds of astrophysical signals / that are sparse in 
some basis, and for which specific prior information is available. 
For each kind of signal, 30 simulations are considered. Observa- 
tions of both kinds of signals are simulated for five hypothetical 
radio interferometers unaffected by instrumental noise, assuming 
that the conditions under which relation ([3]l holds are satisfied. The 
field of view observed on the celestial sphere by the interferometers 
is limited by a Gaussian illumination function A with a full width 
at half maximum (FWHM) of 40 arcminutes of angular opening. 
The original signals considered are defined as sampled images with 
iV = 256 X 256 pixels on a total field of view of 1.8° x 1.8°. 

The first kind of signal consists of a compact object intensity 
field in which the astrophysical objects are represented as a super- 
position of elongated Gaussians of various scales in some arbitrary 
intensity units. The important specific prior information in this case 
is the positivity of the signal. The second kind of signal is of partic- 
ular interest for cosmology. It consists of temperature steps in /^K 
induced by topological defects such as cosmic strings in the zero- 
mean perturbations of the CMB. The string network of interest can 
be mapped as the magnitude of the gradient of the string signal it- 
self. The essential specific prior information in this case resides in 
the fact that the statistical distribution of a string signal may be well 
modelled in wavelet space. One simulation of a compact object in- 
tensity field and the magnitude of the gradient of one simulation of 
a string signal are represented in Figure [T] after multiplication by 
the illumination function. 

As discussed already, we assume uniform random selections 
of visibilities. The five interferometers considered identified by an 
index c with 1 ^ c ^ 5 only differ by their Fourier coverage. This 
coverage is defined by the m/2 randomly distributed frequencies 
probed in one half of the Fourier plane, corresponding to m real 
Fourier coefficients as: m/N — 5c/100. For each configuration, 
the general inverse problem is the one posed in (|5} with the sensing 
matrix $ = defined in ^ if one wishes to undo the multi- 
plication by the illumination function and to recover the original 
signal X. The inverse problem jlOl applies with the sensing matrix 
$ = $,j defined in i ll It if one wishes to recover x. 

For each reconstruction, we compare the performance of the 



Basis Pursuit approaches enhanced by the inclusion of specific 
prior signal information in the minimization problem, with both the 
standard BP^ or BP performance, and the CLEAN performance. As 
the signals considered are sparse or compressible in some basis, we 
do not consider any MEM or WIPE reconstruction, which disregard 
the sparsity information. The performance of the algorithms com- 
pared is evaluated through the signal-to-noise ratio (SNR) of the 
reconstruction for the compact object intensity field, and through 
the SNR of the magnitude of the gradient of the reconstruction for 
the string signal. The SNR of a reconstructed signal s relative to an 
original signal s is technically defined as 

SNR(-^) ^_201ogio^^, (25) 

where o-*-*"*' and a^"^ stand for the sampled standard deviations of 
the residual signal s — s and of the original signal s, respectively. 
It is consequently measured in decibels (dB). 

As far as the computation complexity of the algorithms is 
concerned, notice that both CLEAN and the various Basis Pursuit 
algorithms considered share the same scaling with N at each it- 
eration. This scaling is driven by the complexity of the FFT, i.e. 
0(A'^log7V). The number of iterations required by each algorihm 
is therefore critical in a comparison of computation times. 

4.3 Compact object intensity field 

Each simulation of the compact object intensity field consists of 
100 Gaussians with random positions and orientations, random am- 
plitudes in the range [0, 1] in the chosen intensity units, and ran- 
dom but small scales identified by standard deviations along each 
basis direction in the range [1, 4] in number of pixels. Given their 
structure, such signals are probably optimally modelled by sparse 
approximations in some wavelet basis. But as the maximum possi- 
ble incoherence with Fourier space is reached from real space, we 
chose the sparsity or compressibility basis to be the Dirac basis, 
i.e. 4' = 1^1/2 xjvi/2 • For further simplification of the problem we 
consider the inverse problem l IlOl ) with the sensing matrix $ ^ , for 
reconstruction of the original signal x multiplied by the illumina- 
tion function. 

As no noise is considered, a BP problem is considered in a 
standard compressed sensing approach. However, the prior knowl- 
edge of the positivity of the signal also allows one to pose an en- 
hanced BP-F problem as: 

min I |i subject to J/ = $j,;a;' and x' 0. (26) 

Notice that no theoretical recovery result was yet provided for such 
a problem in the described framework of compressed sensing. But 
the performance of this approach for the problem considered is as- 
sessed on the basis of the simulations. The positivity prior is eas- 
ily incorporated into a convex optim ization solver based on prox- 
imal operator theorv (|Moreaulll962 ). The Douglas-Rachford split- 
ting method jCombettes & Pesguet l2008h guarantees that such an 
additional convex constraint is inserted naturally in an efficient it- 
erative procedure finding the global minimum of the BP-l- problem. 
For simplicity, the stopping criterion of the iterative process is here 
set in terms of the number of iterations: 10''. 

The BP-l- reconstruction of the original signal x reported in 
Figure [T] is also represented in the figure for the configuration 
c — 2. For comparison, the dirty image x'^'''' used in CLEAN and 
obtained by simple application of the adjoint sensing matrix $J. 
to the observed visibilities is also represented. The mean SNR and 
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Figure 1. Top panels: compact object intensity field in some arbitrary intensity units. The original signal multiplied by the illumination function x is reported 
(left), as well as the dirty image x^'*' (center left) and the BP+ reconstruction of x (center right), for the interferometric configuration c = 2. The graph of the 
mean SNR with Ict error bars over 30 simulations is also reported for the CLEAN, BP, and BP+ reconstructions of 5; as a function of the Fourier coverage 
identifying the interferometric configurations (extreme right). Bottom panels: string signal in the CMB in /.tK. The magnitude of the gradient of the original 
signal X re-multiplied by the illumination function is reported (left), as well as the dirty image x^'^^ (center left) and the SBP^ reconstruction of x re-multiplied 
by the illumination function (center right), for the interferometric configuration c = 2. The graph of the mean SNR with Icr eiTor bars over 30 simulations is 
also reported for the CLEAN reconstruction, and for the BP^ and SBP^ reconstructions re-multiplied by the illumination function, as a function of the Fourier 
coverage identifying the interferometric configurations (extreme right). 



corresponding one standard deviation (Icr) error bars over the 30 
simulations are reported in Figiire[T]for the CLEAN reconstruction 
of X with 7 = 0.1, and for the BP and BP-l- reconstructions of x, 
as a function of the Fourier coverage identifying the interferomet- 
ric configurations. All obviously compare very favorably relative to 
the SNR of a;''*', not reported on the graph. One must acknowledge 
the fact that BP and CLEAN provide relatively similar qualities of 
reconstruction. However, the BP reconstruction is actually achieved 
much more rapidly than the CLEAN reconstruction, both in terms 
of number of iterations and computation time. This highlights the 
fact that the BP approach may in general be computationally much 
less expensive. The BP-F reconstruction exhibits a significantly bet- 
ter SNR than the BP and CLEAN reconstructions. The main out- 
come of this analysis thus resides in the fact that the inclusion of the 
positivity prior on the signal significantly improves reconstruction. 
For completeness, let us mention that it was suggested decades ago 
that CLEAN can be unde rstood as some approximatio n of what we 
called the BP-l- approach jMarsh & Richardsonlll987h . 



Notice that the sparsity or compressibility basis is orthonor- 
mal and the error — x*\\2 in the BP reconstruction x* of x is 
theoretically bounded by l |20| ( with e = 0. Assuming saturation of 
this bound, the SNR of the BP reconstruction allows the estimation 
of the maximum sparsity K of the best sparse approximation xk of 
X. Preliminary analysis from the mean SNR of reconstructions over 
the simulations considered suggests that ratios m/ K hold for 
each of the values of m associated with the five interferometric con- 
figurations probed. This result appears to b e in full coherence with 
the accepted empirical ratios quoted above jLustig et al.ll2007t) . 



4.4 String signal in the CMB 

The CMB signal as a whole is a realization of a statistical pro- 
cess. In our setting, the zero-mean temperature perturbations con- 
sidered in /^K may be modelled as a linear superposition of the 
non-Gaussian string signal x made up of steps and of a Gaus- 
sian component g seen as noise. The power spectrum of this as- 
trophysical noise is set by the concordance cosmological model. 
We only include here the so-called primary CMB anisotropics 
jHammond et alJlOOSi) . The typical number, width and spatial dis- 
tribution of long strings or string loops in a given field of view 
are also all governed by the concordance cosmological model. 
Our 30 simulations of the CMB signal are built as a superposi- 
tion of a unique realistic string signal simulation borrowed from 
I 1 ^ ' 

IFraisse et al . (2008) with 30 simulations of the Gaussian correlated 
noise. The string tension p, a dimensionless number related to the 
mass per unit length of string, is up to some extent a free parameter 
of the model. This tension sets the overall amplitude of the signal 
and needs to be evaluated from observations. For the sake of the 
present analysis, we only study the string signal for one realistic 
value p = 3.2 X 10~*, which technically fixes the SNR of the 
observed string signal buried in the astrophysical noise. This value 
is assessed prior to any signal reconstruction, by fitting the power 
spectrum of the data to the sum of th e power spectra of the signal 
and noise on the frequencies probed jHammond et alj |2008). This 
estimation may be considered as very precise at the tension of in- 
terest and is not to be considered as a significant source of error in 
the subsequent reconstruction. 

In this context, preliminary analysis of 1 6 independen t realis - 
tic simulations of a string signal, also from iFraisse et al.l ( l2008h . 
allows one to show that the random process from which the 
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string signal arises is well modelled by GGD's in wavelet space 
jHammond et al.l2008h . We consider a redundant steerable wavelet 
basis with 6 scales j (1 ^ j ^ 6) including low pass and 
high pass axisymmetric filters, and four intermediate scales defin- 
ing steerable wavelets with 6 basis orientations g (1 ^ g ^ 6) 
jSimo ncelli & Freernan 1995). By statistical isotropy, the GGD pri- 
ors TTj for a wavelet coefficient a'^ only depend on the scale: 



TTj (ofuj) oc exp 







Vj- 









(27) 



where w is to be thought of as a multi-index identifying a coeffi- 
cient at given scale j, position i, and orientation q. Assuming in- 
dependence of the wavelet coefficients, the total prior probability 
distribution of the signal is simply the product of the probability 
distributions for each value of w, which reads as 



TT (a) cx exp — 
for a "s" norm 



pUj 



(28) 



(29) 



The exponent parameters Vj are called GGD shape parameters and 
can be considered as a measure of the compressibility of the under- 
lying distribution. Values close to yield very peaked probability 
distributions with heavy tails relative to Gaussian distributions, i.e. 
very compressible distributions. The list of these values at all scales 
reads as: {ui = 0.43, U2 = 0.39, t^s = 0.47, ?;4 = 0.58, us = 
0.76, v% = 1.86}. The signal is thus understood as being well mod- 
elled by a very compressible expansion in its wavelet representation 
and we choose the corresponding redundant basis as the sparsity or 
compressibility basis for the inverse problem: = 'i!^. The list 
values of the GGD scale parameters Uj identifying the variances of 
the distributions at all scales reads as: {mi = 8.9 x 10"'^, U2 = 
2.8 X 10'^U3 = 2.2 X 10^^-U4 = 0.15, W5 = 0.95, ue = 57}. 
In full generality we consider the general inverse problem (O with 
the sensing matrix <1> j , for reconstruction of the original signal x 
non-multiplied by the illumination function. 

Even in the absence of instrumental noise the measured visi- 
bilities thus follow from J16b with a noise term 



n = $ . 



(30) 



representing values of the Fourier transform of the astrophysical 
noise g multiplied by the illumination function. Discarding the very 
local correlations in the Fourier plane introduced by the illumina- 
tion function, one may consider that the measurements are inde- 
pendent and affected by independent Gaussian noise realizations. 
The corresponding noise variance a^. on yr with 1 ^ r ^ m, is 
thus identified from the values of the known power spectrum of g. 

A whitening matrix W^^^^ G R™""" = {(Vi^Jrr' = 
(j~^(5rr'}i^r,r'^m is introduced on the measured visibilities y, so 
that the corresponding visibilities y — Vl^^i, y are affected by in- 
dependent and identically distributed noise, as required to pose a 
BPe problem. This operation corresponds to a matched filtering in 
the absence of which any hope of good reconstruction is vain. A 
BPe problem is thus considered after estimation of p. However, the 
prior statistical knowledge on the signal also allows one to pose 
an enhanced Statistical Basis Pursuit denoise (SBPe) problem. It is 
defined as the minimization of the negative logarithm of the spe- 
cific prior on the signal, i.e. the s norm of the vector of its wavelet 
coefficients, under the measurement constraint: 



Notice that the s norm is similar but still more general than a sin- 
gle Ip norm and no theoretical recovery result was yet provided for 
such a problem in the framework of compressed sensing. Again, 
the performance of this approach for the problem considered is as- 
sessed on the basis of the simulations. Most shape parameters Vj are 
smaller than 1, which implies that the norm defined is not convex. 
We thus reconstruct the signal through the re-we ighted £i norm 
minimization described above /Candes et al."2008'). In this regard, 
we use the SPGLl toolbox ( van den Berg & Friedlander 20083 
The value of in the BP^ and SBP^ problems is taken to be around 

th 9 
the 99 percentile of the x with m degrees of freedom govern- 
ing the noise level estimator. This value also serves as the stopping 
criterion for the CLEAN reconstruction. 

The magnitude of the gradient of the SBPj reconstruction of 
the original signal x reported in Figure[T]is also represented in the 
figure for the configuration c = 2, after re-multiplication by the 
illumination function which sets the field of view of interest. For 
comparison, the magnitude of the gradient of the dirty image a:''*' 
used in CLEAN and obtained by simple application of the adjoint 
sensing matrix ^\ to the observed visibilities is also represented. 

The mean SNR and corresponding one standard deviation 
(la) error bars over the 30 simulations are reported in Figure [T] 
for the CLEAN reconstruction with 7 = 0.1, and for the BP^ and 
SBPe reconstructions re-multiplied by the illumination function, as 
a function of the Fourier coverage identifying the interferometric 
configurations. All obviously compare very favorably relative to 
the SNR of S not reported on the graph. One must still acknowl- 
edge the fact that BP^ and CLEAN provide relatively similar qual- 
ities of reconstruction. The BPe reconstruction is achieved much 
more rapidly than the CLEAN reconstruction, highlighting the fact 
that the BPe approach may in general be computationally much less 
expensive. The SBPe reconstruction exhibits a significantly better 
SNR than the BP and CLEAN reconstructions. 

Let us acknowledge the fact that the re-weighted £1 norm min- 
imization of the SBPe approach proceeds by successive iterations 
of li norm minimization. This unavoidably significantly increases 
the computation time for reconstruction relative to the single f 1 
norm minimization of the BPe approach. Relying on the idea that 
the coefficients of the low pass filter do not significantly participate 
to the identification of the string network itself, our implementa- 
tion of SBPe does not perform any re-weighting at the scale j = 6, 
where — 1 was thus assumed. This restriction allows one to 
keep SBPe computation times similar to those of CLEAN. Let us 
notice however that an even better SNR is obtained by correct re- 
weighting at j = 6, albeit at the cost of a prohibitive increase in 
computation time. 

The main outcome of the analysis is twofold. Firstly, the pres- 
ence of a whitening operation is essential when correlated noise is 
considered. Secondly, the inclusion of the prior statistical knowl- 
edge on the signal also significantly improves reconstruction. 



5 CONCLUSION 

Compressed sensing offers a new framework for image reconstruc- 
tion in radio interferometry. In this context, the inverse problem for 
image reconstruction from incomplete and noisy Fourier measure- 
ments is regularized by the definition of global minimization prob- 
lems in which a generic sparsity or compressibility prior is explic- 



miri a 



subject to \ \y - W^^^^^^-%a'\\2 ^ e. 



(31) 



http://www.cs.ubc.ca/labs/scl/spgll/ 
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itly imposed. These problems are solved through convex optimiza- 
tion. The versatility of this scheme also allows inclusion of specific 
prior information on the signal under scrutiny in the minimization 
problems. We studied reconstruction performances on simulations 
of an intensity field of compact astrophysical objects and of a sig- 
nal induced by cosmic strings in the CMB temperature field, ob- 
served with very generic interferometric configurations. The BPj 
technique provides similar reconstruction performances as the stan- 
dard matching pursuit algorithm CLEAN. The inclusion of specific 
prior information significantly improves the quality of reconstruc- 
tion. 

Further work by the authors along these lines is in preparation. 
In particular, a more complete analysis is being performed to esti- 
mate the lowest string tension down to which compressed sensing 
imaging techniques can reconstruct a string signal in the CMB, in 
more realistic noise and Fourier coverage conditions. In this case, 
given the compressibility of the magnitude of the gradient of the 
string signal itself, TV norm minimization also represents an inter- 
esting alternative to the SBPe problem proposed here. 
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